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We present a formulation of third-order density-functional perturbation theory which is manifestly 
invariant with respect to unitary transfomations within the occupied-states manifold and is partic- 
ularly suitable for a practical implementation of the so called '2n+l' theorem. Our implementation 
is demonstrated with the calculation of the third-order anharmonic coupling coefficients for some 
1 I high-simmetry phonons in Silicon. 
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Density Functional Perturbation Theory (DFPT) is a powerful tool to determine low- 
order derivatives of the ground-state electronic energy of materials with respect to some 
external parameters. The use of DFPT is twofold. On the one hand, it allows to calculate 
response functions — which are directly accessible to experiments — or other measurable 
properties which can be related to response functions — such as e.g. phonon frequencies in 



the adiabatic approximation. On the other hand, it can be used to calculate properties 
of specific, complex, materials which can viewed as small perturbations with respect to 
other, simpler, systems. Since DFPT was demonstrated to be a computationally viable 
technique [1], many applications have appeared belonging to either categories. The first 
^ ', group includes the calculation of elastic constants [2] , dielectric and piezoelectric constants 
[3], and various lattice-dynamical properties [4,5]. Other applications belonging to the 
second group are based on the so called computational alchemy approach to semiconductor 
alloys [6] and superlattices [7] in which the disordered semiconductor is viewed as a small 
perturbation with respect to a reference, periodic system (the virtual crystal) . 

All the applications of DFPT appeared so far are limited to second order in the en- 
ergy. It is a well known result of elementary quantum-mechanics that the knowledge of 
the wavefunction response of a system up to n-th order in the strength of an external 
perturbation is sufficient to determine the energy derivative up to order 2n + 1 [8]. The 
validity of this '2n + 1 theorem' within self-consistent field (SCF) theories has been known 
since several years in the quantum-chemistry community [9] , and recently it has been gen- 
eralized to density-functional theory (DFT) by Gonze and Vigneron [10]. A first important 
conclusion we can draw from this 'theorem' is that the knowledge of the linear response 
of a system to an external perturbation allows to determine the third derivatives of the 
energy with respect to the strength of the perturbation and it gives therefore a practical 
way to link linear and quadratic generalized susceptibilities. The interest in doing so is 
evident: one can in principle obtain higher-order susceptibilities or gain in the accuracy 
achieved by perturbation theory essentially for free. In the following, we will concentrate 
on the formulation by Gonze and Vigneron and will restate it in a form which is free from 
some of its original drawbacks, and is well suited for practical implementations. As an 
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example, we calculate the third-order anharmonic coupling coefficients in Silicon at some 
high-symmetry points of the Brillouin zone (BZ), and compare them with results obtained 
by the frozen-phonon method. 

Let us suppose that the external potential acting on the electrons of a given system 
depends linearly on some external parameter, A: 

v ext (r) = v° ext (r) + Xvl xt (r). (1) 

The case where the dependence of v^ xt upon A is non linear requires a straightforward 
generalization of the results obtained in the linear case and will be considered later. Ac- 
cording to Gonze and Vigneron [10], the third-order derivative of the DFT ground-state 
energy with respect to A reads: 

= 65>Mcf-4IV^>+ / K%v,v\v")n\v)n\v>)n\v")dvdv>dv", (2) 



where the sum runs over occupied {valence) Kohn-Sham (KS) orbitals, ifj 1 and e 1 indi- 
cate the first derivatives of the KS orbitals and energy levels respectively, n 1 and Hg CF 
indicate the corresponding linear corrections to the electron ground-state density and 
KS one-electron hamiltonian, and K 3 is finally the third-order functional derivative of 
the exchange-correlation energy with respect to the electron density: X 3 (r, r',r") = 

. Eq. (2) clearly shows that the calculation of the third-order correc- 
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tion to the energy requires only the knowledge of such ingredients as ifj 1 and e 1 which are 
directly accessible to first-order perturbation theory. 

All the results of DFT must be invariant with respect to unitary transformations of 
the orbitals which do not mix the manifolds of occupied and empty (conduction) states. 
As it stands, Eq. (2) does not manifestly display this invariance. Furthermore, its imple- 
mentation would require the knowledge of the components of the perturbed wavefunction, 
ipl, along all the valence wavefunctions different from ip® itself: {ip^,\ipl) v '^ v . Once again, 
this is innatural because in DFT the variation of any physical property must only depend 
on the variation of the one-electron density matrix which is not affected by components 
of the perturbed valence orbitals along the unperturbed valence manifold. This situation 
is particularly unpleasant when, due to the degeneracy or quasi-degeneracy of some un- 
perturbed valence states, the actual implementation of Eq. (2) would require the use of 
degenerate-state perturbation theory. We will now show that Eq. (2) can be recast in a 
form which requires only the knowledge of the conduction-manifold projection of the ip^s, 
which is manifestly invariant with respect to unitary transformations within the valence 
manifold, and which can be straightforwardly and efficiently implemented using standard 
non-degenerate first-order perturbation theory. 

The second term on the right-hand side (rhs) of Eq. (2) already displays the desired 
unitary invariance, and we concentrate on the first term. Our final result is: 
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where P c = J2 C \^){^\ is the projector over the unperturbed conduction-state manifold 
(from now on, 'c' will indicate an index running over conduction states, while 't> ' indicates 
sums over valence states). Before demonstrating Eq. (3) we notice that it is manifestly 
invariant with respect to unitary transformations within the valence manifold. In fact, 
it is the sum of the trace of a matrix defined over that manifold (first term on the rhs) 
and of the product of two such matrices (second term). The desired invariance derives 
from the invariance property of the trace. The demonstration of Eq. (3) is tedious, but 
straightforward. Let us start from the definition of the first-order correction to the i>-th 
unperturbed valence state, and consider its projections over the valence- and conduction- 
state manifolds: 

\1>l) = Pc\1>l) + Pv\1>l), (4a) 
PM)=J2\^) ^^ , (4b) 

c v c 

v'jiv v v ' 

where P v = 1 — P c is the projector over the valence manifold. Substituting Eq. (4a) into 
the left-hand size of Eq. (3), one obtains the sum of four terms, which will be denoted 
by cc, cv, vc, and vv, according to the couple of projectors appearing inside the matrix 
elements. Inserting Eq. (4c) into the expression of the vv term and separating out terms 
with v' = v" from those with v' ^ v", one obtains: 

V/^ip ^p u/A - V (VWcfWK^ -cl)(<>\Hl CF m 

+ &kA^ (5) 

v^v'^v" ~ e v')( e v - e v") 

Both the terms on the rhs of Eq. (5) vanish because the parities of the numerators and 
those of the denominators with respect to the exchanges v ^ v' v ^ v" are different. 
Let us come now to the cv and vc terms. Using Eqs. (4b) and (4c) and a few algebraic 
manipulations, one obtains: 

Y,{^1\ P o{HIcf - el)Pv\1>l) + (ri\P v {Hl CF - el)PM)) = 

V 

wm CF \ww\Hj, CF \w)(w\Hi, OF M 
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The cc term reads: 

Y.^ 1 v\ P ^ H SCF-tl)PM) = 
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(7) 



The first term on the rhs of Eq. (7) coincides with the first term on the rhs of Eq. (3). 
The second term has the same form as the rhs of Eq. (6), just providing the v = v' terms 
which were missing therein. By combining these terms, we finally obtain: 



V V 

M\ H SCf\'Pc) M H SCfW')('Pv'\ H SCf\'I'v) 



-E 



c,v 



- e° c )(e° v> - eg) 



(8) 



By using Eq. (4b) and the condition that different conduction state are orthogonal to each 
other, we finally arrive at Eq. (3). 

Eq. (3) can be easily generalized to the case where the perturbation depends nonlin- 
early on more than one parameter (as it is actually the case, e.g., in lattice dynamics if A is 
identified with a nuclear displacement) . Suppose there are three different such parameters, 
{Ai, A 2 , A 3 }. Following the notation and the line of reasoning of Ref. [10], we easily arrive 
at the final result: 
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£j\iX2^a _|_ ^A2AiA3 _|_ ^>AiA3A2 _|_ JJJX3X1X2 _|_ JJJX3X2X1 _|_ ^>A2A3Ai ^g-j 
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where 

JJJX1X2X3 _ 

V 

-^(^\Pc\^?)M\H^cfM + 1 J ^ 3 (r,r',r")n Al (r)n A2 (r')n A 3( r ")rfrdr'dr", (10) 

vv' 

where the Aj superscript indicates the derivative with respect to Aj. One sees that when 
the external potential depends lineary on just one parameter, A, Eq. (2) is recovered. 
In the general case where the positions of the nuclei also depend on the A's, one must 
of course add to Eq. (9) the derivative of the ionic contribution to the energy which is 
usually expressed as an Ewald sum. All the ingredients necessary to implement Eq. (10) 
are naturally provided by any computer code aimed at standard second-order DFPT, such 
as the one we routinely use for lattice-dynamical calculations. In the following, we present 
some tests of the above formulation which we have made on the anharmonic coupling 
between lattice distortions of Silicon at selected high-simmetry points of the BZ. 

The equilibrium and lattice-dynamical properties of Silicon have been calculated 
within the local-density approximation, using the plane-wave pseudopotential method. 
We have used the same pseudopotential as in Ref. [4] , plane waves up to a kinetic-energy 
cutoff of 14 Ry, and the (444) Monkhorst-Pack mesh for BZ integrations [11]. Calculations 
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have been done at the F and X points of the BZ both within DFPT and, for comparison, 
by the frozen-phonon method. In the latter case, the unit cell has a lower (rotational 
and/or translational) symmetry, and the set of k-point used for sampling the BZ has been 
modified accordingly. We stress that, as it is the case for the harmonic dynamical matrix 
[4], the calculation of anharmonic coefficients at arbitrary points of the BZ within DFPT 
does not require the use of any supercells, but it only uses wavefunctions and band energies 
calculated for the unperturbed system. There are four independent parameters describ- 
ing the harmonic properties of the crystal within the set of distortions corresponding to 
T and X phonons (the Flto, Xlao, Xta, and Xto frequencies), whereas there are six 
anharmonic constants: one describing the coupling between three T-like phonons, and five 
describing the coupling between one F- and two AT-like phonons. We refer to Ref. [12] for 
a full group-theoretical analysis of the independent coupling coefficients and for an expla- 
nation of the notations we borrow from it. In Table I we compare the third-order coupling 
coefficients calculated in the present work with DFPT and the frozen-phonon method. The 
values obtained with the latter method have been obtained using a procedure analogous to 
the one used in Ref. [12]. As one can see, DFPT give results which are in perfect agreement 
with those obtained by the frozen-phonon method. Actually, they are in principle more 
accurate because DFPT directly provides the energy derivatives without the need of any 
numerical differentiations. 

TABLE I. Comparison of the third-order anharmonic coupling constants between 
phonons at the F and X points of the Brillouin zone in Silicon, as obtained 
by density-functional perturbation theory (DFPT) and the frozen-phonon (FP) 

method. The notations are the same as in Ref. [12]. Units are eV/A 3 . 





B X yz 


I zaa 


-^266 


I zee 


I xac 


I-r- 

ybc 


DFPT 
FP 


295.06 
295.27 


232.41 
232.11 


-35.27 
-35.23 


55.92 
55.44 


447.64 
447.19 


-64.74 
-64.84 



We conclude that DFPT provides an accurate and computationally convenient tool 
for calculating the anharmonic coupling of phonons at arbitrary points of the BZ, with a 
numerical effort which essentially does not depend on the position in the BZ. This opens 
the way to a systematic investigation of such effects in real materials. A calculation of 
the anharmonic-decay phonon lifetimes in semiconductors along the lines presented in this 
paper is presently under way [13]. 

We are grateful to S. de Gironcoli, E. Molinari, and R. Resta for frequent and useful 
discussions. 
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